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The anchor of most integral membrane proteins consists of one or several helices spanning the lipid bilayer. 
The WALP peptide, GWW(LA)^(L)WWA, is a common model helix to study the fundamentals of protein 
insertion and folding, as well as helix-helix association in the membrane. Its structural properties have been 
illuminated in a large number of experimental and simulation studies. In this combined coarse-grained and 
atomistic simulation study, we probe the thermodynamics of a single WALP peptide, focusing on both the 
insertion across the water-membrane interface, as well as folding in both water and a membrane. The potential 
of mean force characterizing the peptide’s insertion into the membrane shows qualitatively similar behavior 
across peptides and three force fields. However, the Martini force field exhibits a pronounced secondary 
minimum for an adsorbed interfacial state, which may even become the global minimum—in contrast to both 
atomistic simulations and the alternative PLUM force field. Even though the two coarse-grained models 
reproduce the free energy of insertion of individual amino acids side chains, they both underestimate its 
corresponding value for the full peptide (as compared with atomistic simulations), hinting at cooperative 
physics beyond the residue level. Folding of WALP in the two environments indicates the helix as the most 
stable structure, though with different relative stabilities and chain-length dependence. 


I. INTRODUCTION 

Transmembrane proteins constitute one of the most 
important biological building blocks, enabling commu¬ 
nication of material and information between a cell 
and its environment, or between different intracellular 
comp art ments.f®^ Despite impressive progress in deter¬ 
mining membrane protein structures,^ aided by techno¬ 
logical advances in fields such as electron tomograph}!^ 
and femtosecond crystallography,l3 the number of known 
structures still lags far behind the case of soluble pro¬ 
teins. Unfortunately, in the absence of structures, the 
options for numerical modeling are limited. This is true 
not only because protein structure prediction remains 
a formidable computational challenge, both for equili¬ 
bration and force-field reasons .1^^^ We also face the ad¬ 
ditional predicament that a lipid bilayer and its sur¬ 
roundings constitute a very highly anisotropic environ¬ 
ment, where everything from dielectric constants to lat¬ 
eral stresses varies dramatically on an Angstrom scale, 
pushing both continuum theory and local thermodynam¬ 
ics to their limits. It should hence not come as a surprise 
that even ostensibly basic questions about structure, lo¬ 
cation and interaction of small peptides in bilayers re¬ 
main difficult to answer .1^ 
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The overwhelming majority of integral membrane pro¬ 
teins is anchored into the lipid bilayer by one or sev¬ 
eral transmembrane a-helices, followed to a much smaller 
fraction by pr oteins where a /3-barrel motif takes over 
that role.^^^^ This is rationalized by the hydrophobic 
environment of the lipid tails, which favors protein con¬ 
formations that minim ize the number of broken backbone 
hydrogen bonds.l^^^^ 

In an effort to better understand membrane proteins 
at a biophysical level, a large body of work has focused 
on studying individual model helices. One common ex¬ 
ample is the sequence of WALP peptides, composed of 
alternating alanine and leucine residues and flanked by 
two tryptophans at each terminus. It was designed to 
resemble a transmembrane helix in membrane pro teins, 
while permitting an easy way to change its length.^^^^^^ 
The arrangement of residues is such that WALP 16 corre¬ 
sponds to the sequence GWW(LA) 5 WWA, while longer 
WALP peptides include more LA repeat units (and oc¬ 
casionally an additional alanine between the final leucine 
and the G-terminal tryptophans). 

Various experimental and simulation studies have shed 
light on the stability of WALP as a transmembrane he¬ 
lix. Experimentally, a combination of NMR methods, 
hydrogen/deuterium exchange, and mass spectrometry 
applied to WALP of different chain lengths, as well as 
lipids of different size, have provided important insight 
into the role of hydrophobic mismatch—the difference 
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between the length of a peptide’s hydrop hobic stretch 
and that of the bilayer’s hydrophobic coreP^^^ For in¬ 
stance, a positive mismatch leads to an average tilt angle 
between the peptide and the membrane normal, a quan¬ 
tity that can be determined from both experiments (e.g., 
quadrupolar splittings from solid-state NMRP^ and 
computer simulationsNotably, Monticelli et al. re¬ 
solved an apparent discrepancy between the average tilt 
angle extracted from experiment versus the same observ¬ 
able calculated in molecular dynamics simulations. Using 
a coarse-grained model, and hence being able to access 
much longer time scales, they showed that both experi¬ 
ment and simulation agree, thus highlighting the impor¬ 
tance of sampling the tilt angle over the microsecond time 
scales relevant for NMR experiment s.^^ Atomistic sim¬ 
ulations later confirmed these findings using enhanced- 
sampling methodologies.!^ 

These studies illuminate the thermodynamics of trans¬ 
membrane helices—not only the stability in the mem¬ 
brane, but also the insertion from water. Using an 
atomistic representation for peptides, but an implicit 
water/membrane model, Im and Brooks showed that 
WALP{ 16,19,23}, starting as an initial random coil, 
would spontaneously insert and fold into a bilawFur¬ 
ther, Nymeyer et al^ and Ulmschneider et al^ demon¬ 
strated insertion and folding of WALP16 in an explicit 
DPPC membrane using enhanced-sampling methodolo¬ 
gies and high-temperature simulations, respectively, to 
alleviate the considerable sampling issues. Some of 
us reported similar findings using PLUM, a recently- 
developed CG model,!^ with and without enhanced 
sampling.!^ 

The potential of mean force (PMF) for the insertion 
of WALP across a water/membrane interface provides 
insight into the thermodynamics of insertion: both in 
terms of the free-energy difference between the two envi¬ 
ronments and the possible existence of intermediate bar¬ 
riers. Structurally, WALP is known to form a helix in 
the membrane, but its conformation in water is largely 
unknown, because its many hydrophobic residues render 
it prone to aggregation at experimentally relevant con¬ 
centrations. Insertion simulations, on the other hand, 
typically work with a single peptide (due to sampling lim¬ 
itations). However, their ability to predict WALP struc¬ 
tures in solution is not merely a matter of the required 
computational resources, but also of the model’s ability 
to describe secondary structure changes in the first place. 
For instance. Bond et al used CG simulations to study 
the thermodynamics of insertion of WALP into a DPPC 
bilayer. Their model, a variant of the CG Martini force 
field,required them to constrain the peptide into a helix 
in all environments,!^ which begs the question whether 
a potential folding/unfolding equilibrium contributes to 
the free energy of insertion. One aim of our present study 
is to address this question. 

The following work investigates the link between 
WALP’s structure and its environment. We rely on the 
CG PLUM force fielcP^to efficiently sample the thermo¬ 


dynamics of insertion across the water-membrane inter¬ 
face, without explicit bias on the secondary structure. To 
gauge the robustness of the results, we carry out equiva¬ 
lent simulations using both the CG Martini force fielcP^, 
bearing in mind its secondary-structure constraints, as 
well as atomistic simulations, despite unavoidable chal¬ 
lenges associated with sampling. While the results agree 
in many qualitative features, we find a number of inter¬ 
esting exceptions which we analyze in some detail. In 
addition, the free-energy profile as a function of helicity 
in both the membrane and water environments provide 
insight into the preferred conformations. 


II. SIMULATION MODELS 
A. Coarse-grained simulations: PLUM force field 

The following describes the CG PLUM force field. The 
associated simulation protocol and parameters used in 
this work are described in Appendix [A} 

The PLUM force field is constructed from the cross- 
par ametr ization of implicit-solvent CG peptid^^ and 
Iipi(p2l33] ijiodels, which we summarize in the following. 

The peptide model includes amino-acid specificity and 
can stabilize different secondary structures using a single 
parametrization, i.e., without explicit bias toward one 
particular conformation. Each amino acid is described 
using four beads: one for the side chain and three for 
the backbone, providing enough resolution to describe 
backbone dihedrals. Phenomenological interactions al¬ 
low the model to reproduce basic properties of peptides 
and proteins, such as excluded volume, hydrophobicity, 
and hydrogen bonds. The model was tuned to qualita¬ 
tively reproduce the Ramachandran plot of tripeptides 
and fold a de novo three-helix bundle. Without chang¬ 
ing the force-field parameters, the model can also sta¬ 
bilize different helical peptides and assemble /3-sheet-rich 
oligomers.!^ The CG model has been appli ed to a variety 
of scenarios involving helical peptides aggregation of 
/3-rich peptides,!^ and /d-barrel formation at the interface 
between virus capsid proteins.!^ 

The lipid model maps a l-palmitoyl-2-oleoyl-sn-gly- 
cero-3-phosphocholine (POPC) lipid into 16 beads, using 
8 bead types to distinguish different chemical moieties.!^ 
Interaction potentials were determined from an iterative- 
Boltzmann inversiorP^ of the radial distribution func¬ 
tions, obtained from an all-atom POPC membrane sim¬ 
ulation. Being an implicit solvent model, the absence of 
water is compensated by a phenomenological attractive 
interaction between tail beads. Free lipids self-assemble 
into a bilayer, which then reproduces elastic properties 
(e.g., the bending modulus), the mass density profile, and 
the orientation of intramolecular bonds .!^ Other neutral 
lipids can be constructed from the set of bead types and 
reach satisfying transferability in terms of structure, area 
per lipid, and temperature dependence of the main phase 
transition.!^ 
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While keeping the individual force-field parameters 
fixed, the cross-parameters between the peptide and 
lipid beads were optimized to reproduce atomistic po¬ 
tential of mean force (PMF) curves of the insertion of 
single amino-acid side chains into a DOPC bilayer 
The cross-parametrization was validated by investigat¬ 
ing a number of structural properties specific to mem¬ 
brane peptides, such as tilt angle, hydrophobic mismatch, 
and transient pore formation from the cooperative ac¬ 
tion of antimicrobial peptidesMore recently, the use 
of a Hamiltonian replica exchange algorithm (more be¬ 
low) assisted in folding several peptides inside the mem¬ 
brane: WALP{ 16,19,23}, as well as the 50-residue-long 
major pVIII coat protein (fd coat) of the filamentous fd 
bacteriophage!^ 


B. Coarse-grained simulations: Martini force field 

Though Martini is a commonly used force field for 
the description of peptide-lipid interactions, we highlight 
some of the key differences with PLUM for completeness. 
The simulation details used throughout this work can be 
found in Appendix [B) 

The coarse-grained Martini model maps on average 
four non-hydrogen atoms into one CG bead, and it can 
describe a wide variety of biomolecules, e.g., water, lip ids, 
proteins, carbohydrates, or small moleculesThe 
key idea is to represent characteristic chemical moyeties 
with a limited set of CG bead types, determined from the 
overall charge, hydrogen-bond capability, and water/oil 
partitioning coefficient 

Martini reproduces a number of lipid-membrane char¬ 
acteristics: self-assembly, area per lipid, elastic prop¬ 
erties, as well as a reasonable bilayer stress profile.®! 
A particularly attractive feature of the model is that 
its building-block approach eases the construction of a 
large variety of molecules, in particular many lipidJ^and 
sterols.®! Due to the mapping of 3 — 4 heavy atoms to 
1 bead there can be some ambiguity with regards to the 
optimal mapping of molecular fragments. For POPC, the 
oleyl tail was originally modeled with 5 beads,®! while 
updated model uses 4 beads.®! 

Martini has been extended t o prot eins, focusing mainly 
on peptide-bilayer interactions.®®! The parametrization 
quite accurately captures the free-energy of the insertion 
of single amino acid side chains and reproduces a number 
of structural properties of model transmembrane helices. 
Though Martini tends to map a similar number of beads 
per amino acid as PLUM, the emphasis is different: a sin¬ 
gle bead represents the backbone while several beads con¬ 
stitute each side chain, providing a better description of 
the sterics. The single-backbone bead description neces¬ 
sitates the use of secondary-structure restraints, present 
in the form of torsional parameters specific to different 
folds (e.g.. Of-helix or /3-sheet). As a result, peptides mod¬ 
eled with Martini cannot (un)fold or refold during the 
simulation. 


C. Atomistic simulations 

The simulation protocol of the atomistic simulations is 
detailed in Appendix 0 

III. RESULTS 

A. Insertion thermodynamics 

Fig.[T] compares the thermodynamics of insertion of a 
single WALP{ 16,19,23} peptide into a POPC membrane 
using different force fields. In each case, the potential of 
mean force (PMF) is displayed as a function of the dis¬ 
tance z between the peptide’s center of mass and that of 
the membrane, which we take as a proxy for its midplane. 
We extended calculations up to 2 : = 5 nm to ensure that 
the peptide was entirely out of the membrane. 

In previous studies, whic h investigated the insertion 
of single amino acids,®®! the bilayer nature was ex¬ 
ploited by simultaneously inserting amino acids into both 
leaflets. This not only increased statistics but also mini¬ 
mized conceivable artifacts due to bilayer asymmetry. In 
our case the size of a peptide makes this strategy unfeasi¬ 
ble, raising the question how the PMF is affected by this 
asymmetric insertion, which stresses the two leaflets dif¬ 
ferently. In Appendix we calculate the resulting elastic 
correction and show it to be negligible. 

Figj^ (a) shows the PMF of WALPI6, WALPI9, and 
WALP23 using the PLUM force field. All curves indi¬ 
cate that the peptide prefers the bilayer over the water 
environment—an expected feature given the hydropho- 
bicity of the amino acids, and in line with the results of 
Bond et a/.®! As we increase the peptide’s length, and 
hence the number of hydrophobic amino acids, the free 
energy of the fully inserted state becomes successively 
smaller. At the bilayer midplane {z = 0), each residue 
contributes on average 1.5 kcal/mol to the free energy of 
insertion. For each PMF, we identify three plateaus: (i) 
close to the bilayer midplane (z ^ 0) the protein samples 
transmembrane conformations; {ii) around the bilayer’s 
interfacial region {z ^ 2 nm) the peptide is still helical, 
but oriented parallel to the surface of the membrane; and 
finally (in) the asymptotic region ( 2 : > 4 nm) where the 
peptide has left the membrane and so its free energy no 
longer depends on z. Representative conformations are 
shown in Fig. for WALPI6, illustrating the transition 
from fully transmembrane to interfacial to desorbed. No¬ 
tice in particular the significant membrane deformation 
occurring ai z ^ 3 nm (see Fig. (d)). It occurs be¬ 
cause the peptide’s free energy gain for staying in contact 
with the membrane outweighs the cost of the concomitant 
elastic deformation—at least for some range of 2 :-values. 
Kopelevich recently showed that these deformations lead 
to an underestimation of the free-energy barrier upon in¬ 
sertion, though the overall free-energy difference should 
be accurate.®! 

Compared to the PLUM results, the PMFs computed 
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FIG. 1. Potential of mean force curves as a function of the distance from the bilayer midplane, z: (a) WALP{16,19,23} from 
the PLUM force field; (b) WALP16 in POPC (“W16”), WALP16 in POPC using the updated force helcP^ with four beads 
for the oleoyl chain (“W16-4B”), WALP16 in DMPC (“W16-D”), and WALP23 in POPC (“W23”) using the Martini model 
with standard water. Note that W16 using Martini’s polarizable water yielded virtually identical results (data not shown); (c) 
WALP16 using the all-atom GROMOS force field. The PLUM energies are mapped from the reduced unit 8 — 0.617 kcal/mol. 
(d), (e), and (f): N-C terminus distance (measured from Ca to Ca) for the respective models. The error of the mean is 
displayed. Note the larger N-C distance for Martini in DMPC (e) close to the bilayer midplane. 


with the Martini model (Fig. (b)) have a noticeably 
different shape. Specifically, all curves exhibit a sec¬ 
ondary minimum corresponding to the interfacial state, 
irrespective of whether the standard or the polarizable 
water model is used (the two curves overlap, only one of 
them is shown). While this interfacial state also exists for 
the PLUM model, as Fig. (c) indicates, its impact on 
the PMF appears much stronger in the Martini model. 
In fact, for WALP16 in a POPC membrane the inter¬ 
facial state is even lower in free energy than the com¬ 
pletely inserted transmembrane state, and hence Martini 
makes a qualitatively different prediction from PLUM 
about thermal equilibrium. In agreement with these re¬ 
sults, a spontaneous transition from transmembrane to 
interfacial states was previously observed by Ramadu- 
rai et using unrestrained simulations of WALP16 in 
lipid membranes made of five or six tailbead-long Martini 
lipids—analogous to the current POPC parametrization. 
In fact, these authors only saw transmembrane-WALP16 
spontaneously transition into the interfacial state when 
they used lipids with long chains. While they did not 
measure a PMF, the barrier from transmembrane to in¬ 
terfacial (Fig. (b); ^ 2 kcal/mol) calculated by us in¬ 
deed suggests the possibility to observe such an event 
spontaneously, given reasonably long simulations. 

To test whether this behavior originates from the 


negative hydrophobic mismatch between WALP16 and 
POPC, we conducted two control simulations: first, we 
kept WALP16 but inserted it into a thinner DMPC bi¬ 
layer; and second, we kept the POPC bilayer but used 
the longer WALP23 peptide. The resulting PMFs (Fig.[^ 
(b)) show that in both cases the transmembrane state be¬ 
comes the most favorable one, even though the interfacial 
state continues to produce a very noticeable metastable 
minimum. This mirrors the observation of Bond et aL, 
who studied WALP23 in DPPC using a customized ver¬ 
sion of the Martini model.l^ 

We measured the membrane thickness from the distri¬ 
bution of distances between the phosphate groups and 
the bilayer midplane projected along the membrane nor¬ 
mal. While PLUM and GROMOS yield similar distri¬ 
butions that peak around 1.8 nm, Martini stabilizes a 
thicker membrane with a peak around 2.1 nm (Fig. ||). To 
compare the impact on the alignment of the peptide, we 
probed the distribution of distances between the trypto¬ 
phan side chains and the bilayer midplane, similarly pro¬ 
jected along the membrane normal. Here again, PLUM 
and GROMOS yield distributions that peak around the 
same point, though the atomistic distribution broadens 
at lower distances. Martini, on the other hand, sam¬ 
ples a distribution shifted by ~ 0.1 nm to higher values. 
The differences in offsets between the tryptophan and 
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FIG. 2. Representative conformations of the WALP16 inser¬ 
tion in POPC using the PLUM force held at different dis¬ 
tances from the bilayer midplane, z: (a) z = 0, (b) z = 1 nm, 
(c) z = 2 nm, (d) z = 3 nm, and (e) z = 4 nm. The peptide 
is depicted in orange, where thick and thin ribbons corre¬ 
spond to the helical and coil states, respectively; the lipids 
are color-coded according to their bead type: purple for the 
hydrocarbon chains, light pastel colors for the interfacial and 
head groups. Rendered with VMD.^ 


phosphate distributions indicate that Martini’s thicker 
membrane will result in the tryptophan side chains being 
buried deeper inside the membrane. Such a deeper inser¬ 
tion will result in a larger energetic penalty, as evid enced 
by the PMF curves of individual side chainj^ ^^ l ^^ l (sam¬ 
pled using the OPLS force field,not GROMOS). 

While the current investigation relied on the origi¬ 
nal POPC Martini model made of five beads for the 
oleoyl chain, Wassenaar et al. recently introduced a 
parametrization using only four beads, thereby reduc¬ 
ing slightly the membrane thickness.!^ The PMF corre¬ 
sponding to the updated force field is shown in Fig. 
“W16-4B”. The reduced hydrophobe mismatch between 
the thinner Martini POPC membrane and WALP16 low¬ 
ers the free energy of the transmembrane state, making it 
roughly equal to that of the inter facial state. This change 
goes into the right direction, but it does not eliminate 
the pronounced minimum of the interfacial state, which 
is absent in the atomistic or PLUM data. This suggests 
that hydrophobic mismatch alone is not the sole reason 
for this feature. 

To explore whether the strong hydrophobic mismatch 
of WALP16 in the 5-bead Martini POPC membrane also 
affects the peptide, we monitored the N- to C-terminal al¬ 
pha carbon distance as a function of 2 : for all force fields 



0 0.5 1 1.5 2 2.5 3 


2 [nm] 


FIG. 3. Probability distributions of normal distances between 
the bilayer midplane and (a) the tryptophan side chains, 
PTiip{z), and (b) the lipid phosphate groups, ppii{z), for 
WALP16 in a POPG membrane modeled by PLUM, Martini, 
and GROMOS. 


(Fig. 0 (d, e, f)) indicates a noticeable stretch for the 
Martini peptide in the region 0 < 2 : < 0.7 nm, which co¬ 
incides with the depth at which the peptide mostly sam¬ 
ples a transmembrane helix (data not shown). WALP16 
in POPC shows a gradual decrease of the N-C distance 
from 2.5 to 2.2 nm between 2 : = 0 and 2 : ~ 2 nm, while 
WALP16 in DMPC displays a sudden drop at 2 : ~ 0.7 nm, 
corresponding to the location of the free-energy barrier 
in Fig. [^(b). Though Martini stabilizes a slightly longer 
helix around ^ = 0, compared to the other force fields, 
the apparent stretching indicates a strong driving force 
to better accommodate a short peptide in the bilayer. On 
the other hand, PLUM and the atomistic simulations (de¬ 
scribed in more details below) do not show any particular 
features close to the bilayer midplane. Unfortunately, the 
atomistic N-C-distance data show a lot of scatter, which 
is clearly a sampling issue. 

For WALP16 in DMPC and WALP23 in POPC, the 
Martini model predicts a pronounced free energy barrier 
(?^ 4 kcal/mol) for the transition from the interfacial to 
the transmembrane state. This is large enough to be¬ 
come a problem in unrestrained simulations that aim to 
study insertion: a peptide which enters the membrane 
from the aqueous phase could get trapped in the interfa¬ 
cial state without transitioning into the transmembrane 
state, even though the latter has a free energy that is 
lower by about 7 kcal/mol. Hall et al. have indeed en¬ 
countered this difficulty during a study that aimed to 
quantify the insertion thermodynamics of various WALP 
peptides in different membranes (using an adapted ver¬ 
sion of Martini).!^ They resorted to co-assembling the 
lipid bilayer in the presence of a WALP peptide and do¬ 
ing statistics of the final state thus obtained (inserted 
or interfacially bound). This protocol suggests that sim¬ 
ply beginning with an interfacially bound peptide was 
not an option, for it would rarely if ever proceed to fully 
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insert—a suspicion which the authors explicitly confirm. 

Despite the rather vivid differences in the shape of the 
PMF, PLUM and Martini largely agree on the free en¬ 
ergy of insertion into the transmembrane state (meaning, 
2 : = 0), provided the hydrophobic mismatch is relaxed. 
This is not completely unexpected, for both models re¬ 
produce the PMFs of in sertio n of single amino-acid side- 
chains into a PC bilayer finding is nevertheless 
nontrivial, because the absolute values do not agree with 
the atomistic ones, as we will discuss below. 

Fig. (c) shows the PMF of WALP16 in POPC, us¬ 
ing the atomistic GROMOS force field. F{z) is largely 
downhill. It exhibits a small shoulder at 2 : = 1 nm, but no 
significant barrier. The location of this shoulder is close 
to the point at which the Martini model finally transi¬ 
tions from interfacial to transmembrane, suggesting that 
this might indeed be the physical origin of this feature, 
but the substantial increase in F{z) by about 10 kcal/mol 
between 1.7 nm and 0.7 nm (observed with Martini) is ab¬ 
sent. Hence, the general shape of the PMF as predicted 
by the PLUM model appears closer to the atomistic data. 

Finally, we wish to point out a curious discrepancy 
between both CG models and the atomistic reference: 
in both GG cases the free energy of WALP16 in its 
equilibrium state (about —20 kcal/mol) is only 2/3 of 
the value predicted in the atomistic simulation (about 
—35kcal/mol). This is surprising, because both models 
capture the free energy of insertion of individual amino 
acids, as predicted atomistically. And while especially in 
the atomistic case one should always be wary of sampling 
issues,!^ and bootstrapping tends to underestimate error 
bars, we do not believe that this is the source of the dis¬ 
crepancy, for it would not suffice to explain a shift by 
15 kcal/mol. Hence, it seems that the difference is real 
and has interesting consequences for modeling. Specifi¬ 
cally, it should be clear that the free energy of insertion 
of an (a-helix consisting of N hydrophobic amino acids 
is not simply the sum of the free energy of insertion of 
each individual amino acid, because there are correla¬ 
tion and cooperativity effects. It seems likely that these 
effects depend not just on the physics captured on the 
coarse-grained level but on more local effects, too. If so, 
CG models of peptides will not capture the insertion free 
energy correctly, even if ostensibly parametrized for pre¬ 
cisely that, and the difference might even be model de¬ 
pendent. Given the large amount of research undertaken 
with these models, it would appear crucial to understand 
this issue better. 


B. Folding in the membrane 

Since the PLUM force field was designed to model 
changes in secondary structure, we can probe the free- 
energy landscape of WALP as a function of helicity— 
using appropriate techniques to ensure accurate sam¬ 
pling. Hamiltonian replica exchange molecular dynamics 
(HREMD) simulations inside the membrane combined 


with the weighted histogram analysis method (WHAM; 
Appendix 0 yields the free energy profiles shown in 
Fig. a Unsurprisingly, the helical state corresponds to 
the free-energy minimumWe note a slight increase in 
the free energy when the helicity approaches 1, illus¬ 
trative of some fraying at the ends of the chain. The 
low-helicity states, on the other hand, are highly sup¬ 
pressed, with free-energy differences ranging from 15 to 
30 kcal/mol. Low but non-zero helicity (~ 0.1) is never 
observed due to the secondary-structure prediction algo¬ 
rithm, which relies on the presence of several (~ 4) con¬ 
secutive amino acids with appropriate hydrogen-bonds 
and dihedrals to assign them in a helical state. We ob¬ 
serve a plateau at low helicity (i.e., 0 — 0.3) followed by 
a sharp, apparently-downhill profile to the helical state. 
Overall, we observe a strong chain-length dependence on 
the free-energy profile. If we plot the three curves against 
the number of broken backbone hydrogen bonds (data 
not shown), the three curves agree more closely in the 
vicinity of their minima, because the change in helic¬ 
ity per broken hydrogen bond depends on the peptide’s 
length. 



FIG. 4. Free energy as a function of helicity for 
WALP{ 16,19,23} inserted in the membrane using the CG 
PLUM force field. Solid lines are mere guides to the eye. The 
grey area roughly indicates the amount of thermal energy at 
body temperature. The conformations underneath illustrate 
an unstructured and helical conformations of WALP 16 in the 
membrane. The PLUM energies are mapped from the reduced 
unit 8 = 0.617 kcal/mol. 
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C. Folding in water 


We then repeat the free energy study as a function of 
helicity from the previous section, but now for WALP 
dissolved in pure water. Fig. (a) shows the free-energy 
profile of the WALP{ 16,19,23} peptides simulated using 
the CG PLUM force field. Just like in the membrane, we 
find a strong preference for helical conformations, with 
free-energy differences between coil and helix in the range 
AF 10 — 25 kcal/mol. 

Interestingly the chain length dependence of the free- 
energy is qualitatively different from the membrane case: 
while WALP 16 again exhibits the lowest free energy at 
any value of the helicity, the profiles for WALP 19 and 
WALP23 are remarkably similar. The qualitative differ¬ 
ence between the two environments is noteworthy, since 
hydrogen bonds are the most likely contributors to the 
free-energy difference. The hydrophobicity will also play 
a larger role in an aqueous environment, as compared to 
the membrane. Yet the interaction strength of hydrogen 
bonds in the model does not depend on whether a bead 
is surrounded by water or lipids .1^ The only noticeable 
difference between the formation of a hydrogen bond in 
water and in the membrane results from the change be¬ 
tween an implicit-water to an explicit-membrane environ¬ 
ment, suggesting an entropic contribution of the model 
itself. Interestingly, we also observed a noticeable change 
in the stability of hydrogen bonds when transferring a 
helix from the water to the membrane environment 
Overall, this behavior may point at a complex interplay 
between the enthalpy (i.e., hydrogen-bonds) and entropy 
of helix formation in waterwhile hydrophobic residues 
immersed in a hydrophobic environment provide more 
straightforward behavior. 

The free-energy profiles shown in Figs. and (a) 
are consistent with the insertion process observed in 


Sec. |III A[ the system sampled a majority of helical con¬ 
formations both in the fully transmembrane state and in 
the aqueous region where the peptide has left the mem¬ 
brane (Fig. |^(a) and (e)). 

We aimed at comparing these findings against refer¬ 
ence atomistic simulations. Using metadynamics,!^ we 
computed the equivalent free-energy profile for WALP 16 
in water (Fig. (b)). The profile shows a minimum at 
80% helicity, which roughly corresponds to 10 over the 
12 possible hydrogen bonds in the peptide, (indicative of 
light fraying of the helix). We observe a fairly complex 
profile with multiple minima, all located above the helical 
state. The helix is therefore the most favorable confor¬ 
mation according to these simulations. Compared to the 
CG results, we find a much narrower profile around the 
minimum. This discrepancy may partially be attributed 
to the difference in defining hydrogen bonds between the 
CG and atomistic simulations (see Appendices [A| and [C|) . 
This also impacts the free energy at low helicity: while 
stride’s definition of a hydrogen bond does not allow us 
to observe any low-helicity conformation (value around 
0.1) in the GG profile, the observable in the metadynam¬ 


ics is continuous along the entire range. 

Aside from difficulties to compare the two curves, we 
point at a possible lack of sampling: the complexity of 
the system makes this free-energy profile difficult to ac¬ 
curately estimate using an atomistic model. The three 
curves shown in Fig. (b), representing the profile after 
simulation times t = 125, 150, and 180 ns per replica, 
are illustrative of the convergence of the profile. We thus 
withhold from further interpreting this curve, and only 
conclude that the helix may indeed be a relevant confor¬ 
mation for WALP 16 in solution. 



FIG. 5. Free energy as a function of helicity for WALP in 
water using (a) CG PLUM simulations on WALP{ 16,19,23} 
and (b) atomistic simulations on WALP 16. Solid lines in (a) 
are mere guides to the eye. The three curves in (b) show the 
profile at different simulation times per replica, i.e., t = 125, 
150, and 180 ns. Note the different scales between (a) and 
(b). The grey area roughly indicates the amount of thermal 
energy at body temperature. The CG conformations under¬ 
neath illustrate random coil and helical states of WALP 16 
in water. The PLUM energies are mapped from the reduced 
unit S = 0.617 kcal/mol. 


IV. CONCLUSIONS 

We performed state-of-the-art thermodynamic calcula¬ 
tions on WALP peptides interacting with a model phos¬ 
pholipid membrane using both coarse-grained (GG) and 
atomistic force fields. The potential of mean force (PMF) 
as a function of penetration depth z indicates increasing 
stability as WALP inserts into the membrane. PLUM 
and GROM OS yield qualitatively similar features: an 
almost-downhill process from water to the fully-inserted 
transmembrane state. Martini, on the other hand, pre¬ 
dicts a distinct minimum for the interfacially bound 













state, which goes along with a pronounced free-energy 
barrier for the transition from the interfacial to the trans¬ 
membrane state for all cases we studied. For WALP16 
in 5-bead POPC this interfacial minimum even becomes 
the global one, in contrast to both PLUM and GRO- 
MOS simulations. Similar behavior was reported in a 
previous Martini study of WALP16 as a function of dif¬ 
ferent lipid-tail sizes, where the transmembrane WALP 
helix spontaneously flipped to the interfacial state in the 
presence of the longer lipids .1^ Though the role of neg¬ 
ative hydrophobic mismatch seems to be predominant 
here, linking this behavior to particular aspects of the 
force field remains difficult. 

Strikingly, PLUM and Martini report very similar 
free-energies of insertion at the bilayer midplane for 
WALP23 in POPC. The agreement likely results from the 
two models’ ability to reproduce the insertion of single 
amino acids in the bilaye r, des pite drastically different 
parametrization strategies.^^^EH Qn the other hand, the 
atomistic GROMOS simulations suggest increased sta¬ 
bility of the transmembrane helix in the membrane: the 
atomistic WALP 16 curve corresponds roughly to the in¬ 
sertion of WALP23 in the CG simulations, yielding a 
discrepancy of ^ 15 kcal/mol. The fact that the two CG 
models underestimate the free energy of insertion by the 
same amount hints at missing correlation and cooper- 
ativity effects beyond the parametrization of individual 
amino acids. This poses questions concerning the coarse- 
graining strategies for peptides on which both PLUM 
and Martini rely, most importantly: under what con¬ 
ditions do matched thermodynamics on the amino acid 
level transfer up to the level of a full peptide? And if it 
does not, what are the dominant sources of discrepancy 
and can we correct for them? 

Though Martini enforces secondary structure, making 
it unable to study the impact of the environment on fold¬ 
ing, plum’s parametrization did allow us to probe this 
behavior in both water and the membrane. Folding in the 
membrane is strongly driven toward the helical state. We 
find roughly linear chain-length dependence on the free- 
energy profile, with longer peptides increasingly penaliz¬ 
ing unstructured random coils. 

Folding in water yielded similar behavior as in the 
membrane, though the chain-length dependence seems 
rather different, likely owing to the complex interplay 
between secondary structure—hydrogen bonds—and ter¬ 
tiary structure— hydrophobicity.^ The results were com¬ 
pared with atomistic metadynamics simulations, which 
also indicate the helix as the most favorable conforma¬ 
tion. Differences in the hydrogen-bond definition, as well 
as sampling difficulties of the atomistic model, make it 
hard to draw further conclusions. Nevertheless, both 
models suggest the helix as a reasonable conformation 
for WALP in water. 

Overall, these findings suggest that enforcing the struc¬ 
ture of a helix throughout the insertion process may rea¬ 
sonably describe the relevant conformational ensemble of 
states. In this sense. Martini’s lack of peptide structural 


rearrangement does not strongly impinge on the results 
for WALP. A better understanding of the contribution 
of (un) folding during the insertion process will require 
the study of a different peptide that shows significantly 
different folds in water and the membrane. A systematic 
comparison of such biomolecular processes using very dif¬ 
ferent computational models (e.g., atomistic vs. coarse¬ 
grained or flexible vs. rigid secondary structure) provides 
a better understanding of the impact of their underlying 
assumption to large-scale thermodynamic properties. 
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Appendix A: PLUM simulation details 

plum’s CG units were constructed from a length 
£ = 1 A, an energy £ = /cBTbody ~ 0.617 kcal/mol 
af ^body = 310 K, and a mass A4. The time unit 
r = MjE ^ 0.1 ps does not properly reflect the dy¬ 
namics of the system, due to the reduction of molecular 
friction during coarse-graining.l^ 

We ran all simulations with the ESPResSo molecular 
dynamics package.l^ A Langevin thermostat and mod¬ 
ified Andersen barostafP^ produced an ensemble with 
constant temperature (T = 1.0 lateral tension 

(E = 0), and vertical box height. Easter integration 
of the equations of motion was achieved using a multi¬ 
timestepping algorithm, setting the short and long time 
steps to = 0.01 r and At = 0.04 r, respect ively.^^ A 
288-POPC lipid membrane was used for all simulations. 
Each peptide was modeled without explicit termini. The 
helicity was determined from the stride secondary- 
structure-prediction algorithm.!^ More sim ulatio n and 
system-setup details are described elsewhere. 

Hamiltonian replica exchange molecular dynamics 
(HREMD]P^ provided enhanced sampling of a peptide in 
both water and the membrane. In particular, we tuned 
the strength of the peptide’s hydrogen-bond interaction, 
i.e., the prefactor e of a modified Lennard-Jones poten¬ 
tial with added directionality (see Appendix l^pl The 
strength of the interaction was modulated by a prefac¬ 
tor, A, where A > 0. We ran HREMD simulations at 
prefactor values from A = 0.1 to A = 1.0, spanning an 
appropriate range of conformational space from fully he¬ 
lical to the complete absence of any helical motif. We 
used 10 and 20 replicas for WALP{16,19} and WALP23, 
respectively. Each replica was run for at least 10^ r. 
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To probe insertion thermodynamics, the distance from 
the bilayer midplane to the peptide was measured from 
the 2 ;-coordinate (i.e., along the bilayer normal) of the 
center of mass of the lipid bilayer to the ^-coordinate of 
the center of mass of the peptide. Umbrella samplin^^ 
restrained the sampled conformational space by restrain¬ 
ing the normal distance between the ^-coordinates of 
the membrane and the peptide. A harmonic restraint of 
spring constant k = 2 £/A was applied at 1 A intervals, 
ensuring enough overlap between the different windows. 
In addition, difficulties associated with sampling PMFs 
of a solute in a lipid membran^^ were addressed here by 
coupling the umbrella sampling with HREMD. Each um¬ 
brella restraint was simulated at 4 interaction prefactors 
A G {0.55,0.70,0.85,1.00} to help sample the conforma¬ 
tional flexibility of the peptide. Each replica was run for 
10^ r, providing an aggregate time of 2 x 10^ r for each 
peptide. 

All enhanced-sampling methodologies were unbiased 
using the weighted histogram analysis method (WHAM; 
see also Appendix TO error bars were calculated 

from bootstrappingP^ 


Appendix B: Martini simulation details 

GROMACS v4.d^ was usedfor the Martini simula¬ 
tions. Martini 2.P^and 2.2FSSIS1I rnodels were used for 
the standard and polarizable water, respectively. A 10-fs 
time step was used, updating the neighbor list every 10 
steps. Lennard-Jones interactions were shifted to zero 
between 0.9 and 1.2 nm. Electrostatic interactions were 
truncated after 1.2 nm with a shifted potential from 0 to 
1.2 nm and a dielectric of 15 (2.5 for polarizable water). 
Although not recommended to be used in classical sim¬ 
ulations, truncation can be used for Martini due to how 
the model has been parametrizedThe temperature of 
310 K was maintained using the V-rescale thermostat!^ 
with a 1-ps time constant. Weak semi-isotropic pressure 
coupling was used with the Berendsen barostat (1-ps time 
constant and 3 x 10“^ bar“^ compressibility)!^. Small bi¬ 
layer patches were simulated with 100 POPC lipids per 
leaflet (121 for the DMPC bilayer). We did not make use 
of specific termini groups. 

We calculated the free energy for transferring a single 
WALP from water to the center of a POPC bilayer (and 
DMPC for Martini 2.2P). A cylindrical position restraint 
was applied from the C^ of the center residue of WALP 
and the center of mass of lipids inside a 1.2-nm-radius 
cylinder centered around the peptide, applied along the 
direction normal to the plane of the bilayer. To prevent 
jumps at the cylinder’s interface, between 1.2 and 1.7 nm 
the weights are switched to zero. A force constant of 500 
kJ/mol/nm^ for the harmonic restraint was used and a 
0.1-nm spacing between adjacent umbrella sampling sim¬ 
ulations, from water (5 nm) to the bilayer center (0 nm). 
Each simulation was run for at least 500 ns. Eree en¬ 
ergy profiles were generated using the weighted histogram 


analysis method (WHAM)PS1 implemented in G_WHA]VP^. 
Error bars were estimated using the bootstrap methocP^ 
with 100 bootstraps. 


Appendix C: Atomistic simulation details 

The final WALP 16 structure from the Martini um¬ 
brella sampling simulations was converted back to atom¬ 
istic representation using the BACKWARD^^ method. 
The GROMOS 54a7 force fields was used on WALP, 
GROMOS on the POPC lipids,EIl and SPCP for wa¬ 
ter. We used a 2-fs time step with bonds to hydro¬ 
gens constrained with the LINGS method.l^ The par¬ 
ticle mesh Ewald summation metho d was used for long- 
range electrostatic interactionsLennard-Jones inter¬ 
actions were shifted from 0.9 to 1.0 nm and truncated 
there after. The V-rescale method^ was used for temper¬ 
ature coupling with a reference temperature of 310 K and 
a 0.1 ps time constant. Pressure was maintained semi- 
isotropically at 1 bar using the Berendsen barostat,!^ a 
2.5-ps time constant and 4.5 x 10“^ bar“^ compressibil¬ 
ity. Eor the umbrella sampling, we increased the har¬ 
monic force constant to 3000 kJ/mol/nm^ and ran each 
simulation for 250 ns. 

To compute the free energy of folding in water for 
WALP 16, a combination of the metadynamic s me thod 
and the parallel tempering scheme was used.l^^^ Ten 
replicas were simulated spanning a temperature range of 
290 — 400 K. An exchange success probability of 9% was 
achieved by applying the well-tempered ensemble (WTE) 
approach, which evenly increases the spread of the po¬ 
tential energy distribution across replic as, while preserv¬ 
ing the same ensemble averagesThe PTMetaD- 
WTE calculations were performed and analyzed with the 
PLUMED2 pluginP 

The PTMetaD-WTE simulations were performed for 
a total period of 180 ns/replica. The relative free-energy 
differences between all of the stable minima were moni¬ 
tored starting after 100 ns/replica and were unchanged 
after this point, so the simulation was terminated after an 
additional 80 ns/replica. The collective variables, 5, bi¬ 
ased in the simulations were a pairwise coordination num¬ 
ber comprising all of the alpha-helical hydrogen bonds 
(i, 2 + 4 pairs) and the peptide’s radius of gyration (alpha 
carbons only; for a discussion on selecting the collective 
variables and a comparison of the different metadynam¬ 
ics techniques, see Ref.^^. The hydrogen-bond collective 
variable is formulated in PLUMED as a summation of 
switching functions of the form: 



with ro = 0.25 nm, n = 6, and m = 9 (note there are 12 
possible (a-helical contacts in WALP 16). These numbers 
were scaled by 12 to provide an approximate “fraction 
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of helicity” in the results. The metadynamics parame¬ 
ters were 0.2 and 0.01 nm for the Gaussian widths of 
the helicity and radius of gyration. Gaussians were de¬ 
posited with a frequency of one per 2 ps, and the con¬ 
vergence of the free-energy estimate was controlled by 
using well-tempered metadynamic^ and a bias factor 
of 10 . As in previous work, an initial simulation period 
(10 ns/replica) was used to equilibrate the replicas to 
a variety of unfolded structures and build up the WTE 
energy bias to achieve overlap between the 10 replicas. 
This simulation used a bias factor of 40 and a Gaussian 
width of 450 kJ/mol. The results presented show a ID 
projection of the 2D metadynamics free-energy surface. 


Appendix D: Estimating free energies from HREMD using 
WHAM 


The weighted histogram analysis method (WHAM) 
provides a minimum variance estimator of the density 
of stat es by combining several simulations of the same 
system.lS^J^ The method is most useful when applied to 
a set of simulations that explore different parts of phase 
space, each contributing to the estimation of thermody¬ 
namic properties of the system. The sampling of phase 
space is enhanced by an appropriate choice of Hamilto¬ 
nians or control parameters (e.g., temperature), which 
together help provide a representative sampling of phase 
space for the process of interest. Though originally ap¬ 
plied to simulations at different temperatures,!^ in the 
following we vary the Hamiltonian of the original sys¬ 
tem, 1-L = Ho + V, where V corresponds to a specific 
part of the Hamiltonian, e.g., an interaction potential. 

In this work, we vary the stre ngth of the protein 
hydrogen-bond interaction potentiaP^^^ 


H(r,'dN,'dc) = ehb 



(Dl) 


J COS^'^N cos^'dc, 

10 


|i?n|,|^?c| <90° 
otherwise 


Each simulation k corresponds to the Hamiltonian H/c = 
Ho + AE, where A > 0 . A = 1 thus corresponds to the 
original Hamiltonian, H, while A 7 ^ 1 alters the propen¬ 
sity to form hydrogen bonds. 

Assuming that all simulations were run at the same 
inverse temperature P = (k^T)~^^ the calculation of the 
free energy as a function of parameter Q is provided by 


PF{Q) oc - In 


E 


Nj exp [/3(Ai - \j)V, - fj] 


(D2) 

where ^(*) bins parameter Q in a discrete set, Nj is the 
number of samples of simulation j, fj is the scaled free 
energy of simulation 7 , i and j sum over simulations, and 
s sums over samples.!^ Determi nation of t he set of fj can 
be obtained by different means.!^* ^^ * ^^ * ^^ * 


Appendix E: Elastic energy of area-leaflet asymmetry upon 
insertion 


Gonsider a bilayer patch that has an area Aq at zero 
tension. If we insert an object into the upper leaflet that 
occupies an area a, the resulting compressive stresses will 
drive an expansion of that leaflet, which in turn puts the 
lower leaflet under tension. In equilibrium, the bilayer 
expands to an area A > Aq, in which a net zero tension 
arises as a balance of compressive and tensile stresses in 
the upper and lower leaflet, respectively. The resulting 
elastic energy contributes to the free energy of insertion 
of the object. How large is it? 

If KA^m = F.aI‘^ is the monolayer stretching modulus, 
the total elastic energy can be written as 


1 {A-A^-af 

^e\ — 2 


. (El) 


Ao 2 ^0 

The still vanishing stress is given by 

'A-Ao-a A-Ao' 




^0 


^0 


(E2) 


from which we find A = Aq + a/ 2 , showing that the area 
mismatch is shared evenly between the two leaflets. The 
total elastic energy is therefore 

Eor WALP, we estimate the area of the inserted object 
as a = 7r(d/2)^, where d = 12 A is the diameter of 
an Qf-helix. In our simulations, we use a relaxed mem¬ 
brane area Aq = (100 A)^, such that (a/Ao)^ ^ 10 “^. 
Given a typical value Ka ~ 250mN/m for the stretch¬ 
ing modulus,!^ we obtain ^ 0.1 /cbT ^ 0.06kcal/mol, 
which is a negligible contribution to the overall free en¬ 
ergy of insertion. 
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